RESULTS NOTIME excluding/regrouping rare species
Running models with different datasets:
all: original data, all species all individuals
exclude: exclude rare species
regroup: regrouping observations of rare species in ONE species group
Using only the 5x5 m quadrat scale
#Figures colors/labels
cores <- c("#c00000", "#b666d2", "#70ad47", "#ffc000")
labvit <- as_labeller(c(grow = "Growth", mort="Mortality",
rec="Recruitment"))
grplab = as_labeller(c(quadrat = "Space", `quadrat:sp`= "Space X Species",
sp="Species", Residual = "residual"))
grpvit <- as_labeller(c(grow = "Growth", mort="Mortality",
rec="Recruitment",
quadrat = "Space", `quadrat:sp`= "Species X Space",
sp="Species", Residual = "Residual"))Data
Growth
grp <- c("all","exclude", "regroup")
tudo <- list()
path1 = here("models_outputs/models_notime/")
for (i in 1:length(grp)){
path <- paste0(path1, grp[i], "/grow/table")
files = list.files(path)
all <- lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
names(all) <- substr(files,1,nchar(files)-6)
tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")
#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$fplot,tes$data)/5##
## all exclude regroup
## ama 1 1 1
## bci 6 6 6
## edo 2 2 2
## fus 3 3 3
## idc 1 1 1
## kor 1 1 1
## lam 3 3 3
## ldw 1 1 1
## len 2 2 2
## lpl 1 1 1
## luq 3 3 3
## mos 2 2 2
## pas 4 4 4
## scbi 1 1 1
## serc 1 1 1
## sin 2 2 2
## ucsc 2 2 2
## wab 2 2 2
## wfdp 1 1 1
## wyw 3 3 3
## zof 1 1 1
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"
grow <- tes %>% filter(term != "intercept") %>%
left_join(intercepts[,c(1:4,9)], by= c("id", "data", "fplot", "time")) %>%
group_by(id, data, fplot, time) %>%
mutate(VPC = variance/sum(variance)) %>% ungroup() %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual")) %>%
as.data.frame()
# divide the sd of the terms by the mean growth (intercept) for each dataset.
grow$stdev = grow$Estimate/grow$interceptMort
grp <- c("all","exclude", "regroup")
tudo <- list()
path1 = here("models_outputs/models_notime/")
for (i in 1:length(grp)){
path <- paste0(path1,grp[i], "/mort/table")
files = list.files(path)
all <- lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
names(all) <- substr(files,1,nchar(files)-6)
tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")
#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$fplot,tes$data)/5##
## all exclude regroup
## ama 1 1 1
## bci 7 7 7
## edo 2 2 2
## fus 3 3 3
## idc 1 1 1
## kor 1 1 1
## lam 3 3 3
## ldw 1 1 1
## len 2 2 2
## lpl 1 1 1
## luq 3 3 3
## mos 2 2 2
## pas 5 5 5
## scbi 1 1 1
## serc 1 1 1
## sin 2 2 2
## ucsc 2 2 2
## wab 2 2 2
## wfdp 1 1 1
## wyw 3 3 3
## zof 1 1 1
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"
mort <- tes %>% filter(term != "intercept") %>%
left_join(intercepts[,c(1:4,9)], by= c("id", "data", "fplot", "time")) %>%
group_by(id, data, fplot, time) %>%
mutate(VPC = variance/sum(variance)) %>% ungroup() %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual")) %>%
as.data.frame()
mort$stdev = mort$EstimateRec
grp <- c("all","exclude", "regroup")
tudo <- list()
path1 = here("models_outputs/models_notime/")
for (i in 1:length(grp)){
path <- paste0(path1,grp[i], "/rec/table")
files = list.files(path)
all <- lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
names(all) <- substr(files,1,nchar(files)-6)
tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")
#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$fplot,tes$data)/5##
## all exclude regroup
## bci 7 7 7
## edo 2 2 2
## fus 3 3 3
## idc 1 1 1
## kor 1 1 1
## lam 3 3 3
## ldw 1 1 1
## len 2 2 2
## lpl 1 1 1
## luq 3 3 3
## mos 2 2 2
## pas 5 5 5
## scbi 1 1 1
## serc 1 1 1
## sin 2 2 2
## ucsc 2 2 2
## wab 2 2 2
## wfdp 1 1 1
## wyw 3 3 3
## zof 1 1 1
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"
rec <- tes %>% filter(term != "intercept") %>%
left_join(intercepts[,c(1:4,9)], by= c("id", "data", "fplot", "time")) %>%
group_by(id, data, fplot, time) %>%
mutate(VPC = variance/sum(variance)) %>% ungroup() %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual")) %>%
as.data.frame()
rec$stdev = rec$EstimateExcluding recruitment for Wytham Woods, as in “all data”.
rec <- rec %>% filter(fplot!="wyw")Combining results.
comb <- bind_rows(list(grow=grow, mort=mort, rec=rec), .id="vital") %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp","Residual")) %>%
ungroup() Comparing groups
comb %>%
ggplot(aes(x=data, y=VPC)) +
geom_col(aes(fill=term, color=term), alpha=0.7) +
scale_fill_manual(values=cores) +
scale_color_manual(values=cores) +
facet_grid(vital~fplot+time) +
theme(axis.text.x = element_text(angle=45, hjust=1))mean over fplots
comb %>% group_by(vital,data,term) %>% summarise(VPC=mean(VPC)) %>%
ggplot(aes(x=data, y=VPC)) +geom_col(aes(fill=term, color=term), alpha=0.7) +
scale_fill_manual(values=cores) +
scale_color_manual(values=cores) +
facet_grid(~vital)+
theme(axis.text.x = element_text(angle=45, hjust=1))+
labs(tag="a)")#ggsave("figs/FIG_S4.4a_groups_mean.jpg", width=20, height = 10, units = "cm")comb %>%
ggplot(aes(x=data,y=VPC, col=data)) + geom_boxplot()+
facet_grid(vital~term) +
theme(axis.text.x = element_text(angle=45,hjust=1),
legend.position = "none")comb %>% mutate(xis = as.numeric(as.factor(data))) %>%
ggplot(aes(x=xis,y=stdev, col=paste(fplot,time))) + geom_line() +
facet_grid(vital~term) +
scale_x_continuous(breaks = 1:3, labels=levels(comb$data),
name="") +
theme(axis.text.x = element_text(angle=45, hjust=1),
legend.position = "none")comb %>% mutate(xis = as.numeric(as.factor(data))) %>%
ggplot(aes(x=xis,y=VPC, col=paste(fplot,time))) + geom_line() +
facet_grid(vital~term) +
scale_x_continuous(breaks = 1:3, labels=levels(comb$data),
name="") +
theme(axis.text.x = element_text(angle=45, hjust=1),
legend.position = "none")Difference to all data
Absolute difference
dif <- comb %>%
pivot_wider(id_cols = c(vital,fplot,time,term),
names_from = data,
values_from = c(stdev,VPC, richness)) %>%
mutate(stdev.dif_exclude = stdev_exclude - stdev_all,
VPC.dif_exclude = VPC_exclude - VPC_all,
stdev.dif_regroup = stdev_regroup - stdev_all,
VPC.dif_regroup = VPC_regroup - VPC_all) %>%
select(vital,fplot,time,term,richness_exclude,
stdev.dif_exclude,
stdev.dif_regroup, VPC.dif_exclude, VPC.dif_regroup) %>%
pivot_longer(6:9, names_to = c("index", "data"), names_sep = "_",
#names_pattern = "(.*_).(.*)",
values_to = "value") %>%
pivot_wider(names_from = index, values_from = value)dif %>%
ggplot(aes(x=data, y=stdev.dif)) +
facet_grid(vital~term, scales= "free_x") +
geom_boxplot() +
geom_quasirandom(aes(col=fplot)) +
ylab("SD no.rare - SD all") + xlab("") +
geom_hline(yintercept = 0, linetype= "dotted") +
ggtitle("Difference in SD") +
theme(legend.position = "right",
panel.background = element_rect(colour="lightgray"),
axis.text.x = element_text(angle=45, hjust=1))dif %>%
ggplot(aes(x=data, y=VPC.dif)) +
facet_grid(vital~term, scales= "free_x") +
geom_boxplot() +
geom_quasirandom(aes(col=fplot)) +
ylab("VPC no.rare - VPC all") + xlab("") +
geom_hline(yintercept = 0, linetype= "dotted") +
ggtitle("Difference in VPC") +
theme(legend.position = "right",
panel.background = element_rect(colour="lightgray"),
axis.text.x = element_text(angle=45, hjust=1))Relative difference SD
difr <- comb %>%
pivot_wider(id_cols = c(vital,fplot,time,term),
names_from = data,
values_from = c(stdev,VPC, richness)) %>%
mutate(stdev.dif_exclude = (stdev_exclude - stdev_all)*100/stdev_all,
stdev.dif_regroup = (stdev_regroup - stdev_all)*100/stdev_all) %>%
select(vital,fplot,time,term,richness_exclude,
stdev.dif_exclude,stdev.dif_regroup) %>%
pivot_longer(6:7, names_to = "data",
values_to = "stdev.dif") %>%
mutate(data = substr(data,11, nchar(data)))difr %>%
ggplot(aes(x=data, y=stdev.dif)) +
facet_grid(vital~term, scales= "free_x") +
geom_boxplot() +
geom_quasirandom(aes(col=fplot)) +
ylab("Proportional difference (%)") + xlab("") +
geom_hline(yintercept = 0, linetype= "dotted") +
ggtitle("Relative difference in SD") +
theme(legend.position = "right",
panel.background = element_rect(colour="lightgray"),
axis.text.x = element_text(angle=45, hjust=1))difr %>% group_by(vital, term, data) %>% summarise(mean=mean(stdev.dif, na.rm=T)) %>%
pivot_wider(names_from = term, values_from = mean)## # A tibble: 6 × 6
## # Groups: vital [3]
## vital data quadrat `quadrat:sp` sp Residual
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 grow exclude 11.0 8.47 -11.8 13.2
## 2 grow regroup 8.71 3.97 -12.7 13.9
## 3 mort exclude 0.0278 -1.03 -14.9 0
## 4 mort regroup -2.62 -0.591 -16.4 0
## 5 rec exclude -0.565 -0.299 -9.26 0
## 6 rec regroup 0.512 -0.768 -10.4 0
table
Relative differences in SD
difr %>% group_by(vital, term, data) %>%
summarise(mean=mean(stdev.dif, na.rm=T)) %>%
pivot_wider(names_from = c(vital,data), values_from = mean) %>%
kable(digits=1)| term | grow_exclude | grow_regroup | mort_exclude | mort_regroup | rec_exclude | rec_regroup |
|---|---|---|---|---|---|---|
| quadrat | 11.0 | 8.7 | 0.0 | -2.6 | -0.6 | 0.5 |
| quadrat:sp | 8.5 | 4.0 | -1.0 | -0.6 | -0.3 | -0.8 |
| sp | -11.8 | -12.7 | -14.9 | -16.4 | -9.3 | -10.4 |
| Residual | 13.2 | 13.9 | 0.0 | 0.0 | 0.0 | 0.0 |
Mean absolute differences in VPC
dif %>% group_by(vital, term, data) %>%
summarise(mean=mean(VPC.dif, na.rm=T)) %>%
pivot_wider(names_from = c(vital,data), values_from = mean) -> write
kable(write)| term | grow_exclude | grow_regroup | mort_exclude | mort_regroup | rec_exclude | rec_regroup |
|---|---|---|---|---|---|---|
| quadrat | 0.0048918 | 0.0039062 | 0.0081885 | 0.0061325 | 0.0069821 | 0.0130636 |
| quadrat:sp | 0.0079957 | 0.0028262 | 0.0139204 | 0.0162427 | 0.0076814 | 0.0044177 |
| sp | -0.0831334 | -0.0908756 | -0.0704275 | -0.0754252 | -0.0321646 | -0.0368164 |
| Residual | 0.0702459 | 0.0841432 | 0.0483186 | 0.0530499 | 0.0175012 | 0.0193351 |
Comparing with NUMBER of rare species
dif %>%
ggplot(aes(x=richness_exclude, y=stdev.dif, shape=data)) +
geom_point(aes(col=fplot)) +
geom_smooth(method="lm", se=F, aes(linetype=data)) +
facet_grid(vital~term) +
ylab("Difference in SD to original data") +
xlab("log N of rare species") +
geom_hline(yintercept = 0, linetype= "dotted") +
ggtitle("Difference in SD") +
scale_x_log10()+
theme(axis.text.x = element_text(angle=45, hjust=1),
)difr %>%
ggplot(aes(x=richness_exclude, y=stdev.dif, shape=data)) +
geom_point(aes(col=fplot)) +
geom_smooth(method="lm", se=F, aes(linetype=data)) +
facet_grid(vital~term) +
ylab("Relative difference in SD to original data") +
xlab("log N of rare species") +
geom_hline(yintercept = 0, linetype= "dotted") +
ggtitle("Proportional difference in SD to original data") +
theme(axis.text.x = element_text(angle=45, hjust=1),
panel.background = element_rect(fill="white", color="black"))pm <- dif %>%
ggplot(aes(x=richness_exclude, y=VPC.dif, shape=data)) +
geom_point(aes(col=fplot)) +
facet_grid(vital~term) +
geom_smooth(method="lm", se=F, aes(linetype=data)) +
ylab(" Absolute difference in VPC to original data") +
xlab("log N of rare species") +
geom_hline(yintercept = 0, linetype= "dotted") +
theme(axis.text.x = element_text(angle=45, hjust=1),
panel.background = element_rect(fill="white", color="black"))
pmDirichlet regression for exclude results with species richness
Richness by rarefaction at the 6ha plot size (sampling unit 20x20m)
load(here("data", "rarefaction.curves.Rdata"))
rich.rare <- rare[rare$sites ==150,2:3]
colnames(rich.rare)[1] <- "richness.rarefaction"
comb <- comb %>% left_join(rich.rare, by="fplot")Latitude
load(here("data", "plots_structure.Rdata"))
comb <- comb %>% left_join(plots.structure[,1:2], "fplot")#mean for each forest
mcomb <- comb %>% group_by(vital,data, fplot,term, lat) %>%
summarise(VPC = mean(VPC),
stdev = mean(stdev),
richness.rarefaction = mean(richness.rarefaction))GROWTH
predirig <- mcomb %>%
filter(vital == "grow", data == "exclude") %>%
select(fplot, term, VPC, richness.rarefaction,lat)
diridatag <- predirig %>%
pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
mutate(log.rich = log(richness.rarefaction),
log.rich.o=log.rich) %>%
mutate_at(vars( log.rich), scale)
vpcg <- DR_data(diridatag[,c("quadrat", "sp", "quadrat:sp", "Residual")])
#plot(vpc)grow2 <- DirichReg(vpcg~log.rich, data=diridatag, model="alternative", base=4)
summary(grow2)## Call:
## DirichReg(formula = vpcg ~ log.rich, data = diridatag, model = "alternative",
## base = 4)
##
## Standardized Residuals:
## Min 1Q Median 3Q Max
## quadrat -0.9058 -0.4992 -0.2701 -0.0166 0.9614
## sp -1.5134 -0.7438 -0.3542 0.2168 3.7107
## quadrat:sp -1.7533 -0.5644 -0.2376 0.3694 2.6473
## Residual -2.0096 -0.8659 0.3287 1.1003 1.8630
##
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6205 0.1744 -15.026 < 2e-16 ***
## log.rich 0.4740 0.1675 2.829 0.00466 **
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.09746 0.09944 -11.036 <2e-16 ***
## log.rich -0.20184 0.10053 -2.008 0.0447 *
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.44533 0.11231 -12.869 <2e-16 ***
## log.rich -0.08056 0.11309 -0.712 0.476
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
##
## PRECISION MODEL:
## ------------------------------------------------------------------
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.4105 0.1789 19.06 <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood: 102.9 on 7 df (38 BFGS + 1 NR Iterations)
## AIC: -191.7, BIC: -184.4
## Number of Observations: 21
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
newdata <- expand.grid(log.rich = seq(min(diridatag$log.rich),
max(diridatag$log.rich), length.out=10))
pred <- predict(grow2, newdata = newdata, se=T)
confint(grow2)##
## 95% Confidence Intervals (original form)
##
## - Beta-Parameters:
## Variable: quadrat
## 2.5% Est. 97.5%
## (Intercept) -2.962 -2.621 -2.279
## log.rich 0.146 0.474 0.802
##
## Variable: sp
## 2.5% Est. 97.5%
## (Intercept) -1.292 -1.097 -0.903
## log.rich -0.399 -0.202 -0.005
##
## Variable: quadrat:sp
## 2.5% Est. 97.5%
## (Intercept) -1.665 -1.445 -1.225
## log.rich -0.302 -0.081 0.141
##
## Variable: Residual
## variable omitted
##
## - Gamma-Parameters
## 2.5% Est. 97.5%
## (Intercept) 3.06 3.41 3.76
colnames(pred) <- colnames(vpcg)
newdata$log.rich.o <- newdata$log.rich*sd(diridatag$log.rich.o) +
mean(diridatag$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredg <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")
newpredg %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
"Residual")) %>%
ggplot(aes(x=rich.o, y=pred, col=term)) +
geom_line() + facet_grid(~term) +
geom_point(data=predirig, aes(x=richness.rarefaction, y=VPC, col=term)) +
scale_x_log10() +
ggtitle("Mod log.rich")Residuals
resid <- residuals(grow2, type = "standardized")
resid <- data.frame(resid[,1:4])
resid$log.rich <- diridatag$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")
ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
facet_grid(~term)MORT
predirim <- mcomb %>%
filter(vital == "mort", data == "exclude") %>%
select(fplot, term, VPC, richness.rarefaction,lat)
diridatam <- predirim %>%
pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
mutate(log.rich = log(richness.rarefaction),
log.rich.o=log.rich) %>%
mutate_at(vars( log.rich), scale)
vpcm <- DR_data(diridatam[,c("quadrat", "sp", "quadrat:sp", "Residual")])
#plot(vpc)mort2 <- DirichReg(vpcm~log.rich,data=diridatam, model="alternative", base=4)
summary(mort2)## Call:
## DirichReg(formula = vpcm ~ log.rich, data = diridatam, model = "alternative",
## base = 4)
##
## Standardized Residuals:
## Min 1Q Median 3Q Max
## quadrat -1.2846 -0.6061 -0.0454 0.4078 1.8812
## sp -1.9682 -0.8456 0.1570 0.9680 2.6456
## quadrat:sp -1.8027 -0.6593 -0.0102 0.2449 1.6414
## Residual -2.9957 -0.2679 0.0598 0.6229 1.6857
##
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.6901 0.1375 -12.295 < 2e-16 ***
## log.rich 0.4103 0.1473 2.786 0.00534 **
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.87149 0.10008 -8.708 <2e-16 ***
## log.rich -0.14106 0.09865 -1.430 0.153
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2178 0.1135 -10.73 <2e-16 ***
## log.rich -0.2471 0.1217 -2.03 0.0424 *
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
##
## PRECISION MODEL:
## ------------------------------------------------------------------
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.3788 0.1756 19.24 <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood: 88.64 on 7 df (43 BFGS + 1 NR Iterations)
## AIC: -163.3, BIC: -156
## Number of Observations: 21
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
newdata <- expand.grid(log.rich = seq(min(diridatam$log.rich),
max(diridatam$log.rich), length.out=10))
pred <- predict(mort2, newdata = newdata, se=T)
confint(mort2)##
## 95% Confidence Intervals (original form)
##
## - Beta-Parameters:
## Variable: quadrat
## 2.5% Est. 97.5%
## (Intercept) -1.960 -1.69 -1.421
## log.rich 0.122 0.41 0.699
##
## Variable: sp
## 2.5% Est. 97.5%
## (Intercept) -1.068 -0.871 -0.675
## log.rich -0.334 -0.141 0.052
##
## Variable: quadrat:sp
## 2.5% Est. 97.5%
## (Intercept) -1.440 -1.218 -0.995
## log.rich -0.486 -0.247 -0.009
##
## Variable: Residual
## variable omitted
##
## - Gamma-Parameters
## 2.5% Est. 97.5%
## (Intercept) 3.04 3.38 3.72
colnames(pred) <- colnames(vpcm)
newdata$log.rich.o <- newdata$log.rich*sd(diridatam$log.rich.o) +
mean(diridatam$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredm <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")
newpredm %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
"Residual")) %>%
ggplot(aes(x=rich.o, y=pred, col=term)) +
geom_line() + facet_grid(~term) +
geom_point(data=predirim, aes(x=richness.rarefaction, y=VPC, col=term)) +
scale_x_log10() +
ggtitle("Mod log.rich")Residuals
resid <- residuals(mort2, type = "standardized")
resid <- data.frame(resid[,1:4])
resid$log.rich <- diridatam$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")
ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
facet_grid(~term)REC
predirir <- mcomb %>%
filter(vital == "rec", data == "exclude") %>%
select(fplot, term, VPC, richness.rarefaction,lat)
diridatar <- predirir %>%
pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
mutate(log.rich = log(richness.rarefaction),
log.rich.o=log.rich) %>%
mutate_at(vars( log.rich), scale)
vpcr <- DR_data(diridatar[,c("quadrat", "sp", "quadrat:sp", "Residual")])
#plot(vpc)rec2 <- DirichReg(vpcr~log.rich,data=diridatar, model="alternative", base=4)
summary(rec2)## Call:
## DirichReg(formula = vpcr ~ log.rich, data = diridatar, model = "alternative",
## base = 4)
##
## Standardized Residuals:
## Min 1Q Median 3Q Max
## quadrat -1.5124 -0.6320 -0.2358 0.8198 2.1508
## sp -2.8666 -0.2882 0.1438 0.7867 1.6337
## quadrat:sp -1.4254 -0.6122 -0.2100 0.3186 2.8831
## Residual -1.5080 -0.6407 -0.0523 0.5100 3.1244
##
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.586194 0.129600 -4.523 6.09e-06 ***
## log.rich 0.005999 0.140757 0.043 0.966
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.08768 0.11441 -0.766 0.443
## log.rich -0.71600 0.11995 -5.969 2.38e-09 ***
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.70344 0.13331 -5.277 1.32e-07 ***
## log.rich -0.07571 0.14027 -0.540 0.589
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
##
## PRECISION MODEL:
## ------------------------------------------------------------------
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.2418 0.1825 17.77 <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood: 71.27 on 7 df (43 BFGS + 1 NR Iterations)
## AIC: -128.5, BIC: -121.9
## Number of Observations: 19
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
newdata <- expand.grid(log.rich = seq(min(diridatar$log.rich),
max(diridatar$log.rich), length.out=10))
pred <- predict(rec2, newdata = newdata, se=T)
colnames(pred) <- colnames(vpcr)
newdata$log.rich.o <- newdata$log.rich*sd(diridatar$log.rich.o) +
mean(diridatar$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredr <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")
newpredr %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
"Residual")) %>%
ggplot(aes(x=rich.o, y=pred, col=term)) +
geom_line() + facet_grid(~term) +
geom_point(data=predirir, aes(x=richness.rarefaction, y=VPC, col=term)) +
scale_x_log10() +
ggtitle("Mod log.rich")Residuals
resid <- residuals(rec2, type = "standardized")
resid <- data.frame(resid[,1:4])
resid$log.rich <- diridatar$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")
ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
facet_grid(~term)FIGURE
Calculated prediction interval
source(here("scripts/prediction_intervals_dirichlet_exclude.R"), local=T)
#load("results/prediction_intervals_dirichlet_exclude.Rdata") # quantspredis <- bind_rows(list(grow=newpredg, mort=newpredm, rec=newpredr),
.id="vital") %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp","Residual"))
pontos <- bind_rows(list(grow=predirig, mort= predirim, rec=predirir), .id="vital")%>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))
quants <- bind_rows(list(grow=quantg,mort=quantm, rec=quantr), .id="vital")%>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))# pvalues
res <- data.frame(vital = c("grow", "mort", "rec"),
P = NA,
term = "Residual")
test <- summary(grow2)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsg <- vals[ grep("log.rich", vals$coef),] %>%
mutate(term=c("quadrat", "sp", "quadrat:sp"))
test <- summary(mort2)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsm <- vals[ grep("log.rich", vals$coef),] %>%
mutate(term=c("quadrat", "sp", "quadrat:sp"))
test <- summary(rec2)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsr <- vals[ grep("log.rich", vals$coef),] %>%
mutate(term=c("quadrat", "sp", "quadrat:sp"))
pvals <- bind_rows(list(grow=valsg,mort=valsm,rec=valsr), .id="vital") %>%
dplyr::select(vital,`Pr(>|z|)`, term) %>% rename(P = `Pr(>|z|)`)
pvals <- rbind(pvals,res) %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))
pvals$x <- 300
pvals$y <- 0.70
pvals$sig <- paste0("p = ", round(pvals$P,3))
#pvals$sig[pvals$P >0.0166] <- "ns"
pvals$sig[pvals$term=="Residual"] <- ""
pvals$sig[pvals$sig=="p = 0"] <- "p < 0.0001"library(wesanderson) #better colour palette
# Gradient color para latitutde
pal <- wes_palette("Zissou1",20, type = "continuous")[20:1] # azul frio-temperadofdiri_lat <-ggplot(predis, aes(x=rich.o, y=pred)) +
geom_line(size=1) +
facet_grid(vital~term, labeller=grpvit) +
geom_smooth(data=quants,aes(x=rich.o, y=lower), se=F, size=0.1)+
geom_smooth(data=quants,aes(x=rich.o, y=upper), se=F, size=0.1) +
geom_ribbon(data=quants,aes(x=rich.o, ymin=lower, ymax=upper,
y=mean), alpha=0.05, fill="blue",
size=0)+
geom_point(data=pontos, aes(x=richness.rarefaction, y=VPC, fill=abs(lat)),
col="black", size=2.5, pch=21)+
scale_fill_gradientn(colors=pal) +
scale_y_continuous(minor_breaks = NULL) +
scale_x_log10() +
geom_text(data=pvals, aes(x=x,y=y, label=sig), size=3, hjust=1,vjust=0,
col="black")+
xlab("Species richness (log10 scale)")+
ylab("VPC") +
theme_cowplot() +
theme(panel.background = element_rect(colour="lightgray"),
legend.position = "none")
fdiri_latpng(here("figures/FIG_S5.5_dirichlet_regressions_lat_exclude.png"), height = 600, width=800, res=100)
fdiri_lat
dev.off()## quartz_off_screen
## 2